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ABSTRACT 

o 

■ We have developed a new massively-parallel radiation-hydrodynamics code (Cosmos) for New- 

toman and relativistic astrophysical problems that also includes radiative cooling, self-gravity, 
and non-equilibrium, multi-species chemistry. Several numerical methods are implemented for 
the hydrodynamics, including options for both internal and total energy conserving schemes. Ra- 
diation is treated using flux-limited diffusion. The chemistry incorporates 27 reactions, including 
both collisional and radiative processes for atomic hydrogen and helium gases, and molecular 
hydrogen chains. In this paper we discuss the equations and present results from test problems 
carried out to verify the robustness and accuracy of our code in the Newtonian regime. An earlier 
Qs, ' paper presented tests of the relativistic capabilities of Cosmos. 

O ' 

\ Subject headings: diffusion — hydrodynamics — instabilities — methods: numerical — shock 

o 
m 
o 

' 1. Introduction 

Oh; 

This is the second in a series of papers using a new code (Cosmos) that we have developed for a 
> ■ broad range of astrophysical problems, including scalar- and radiation- field dominated processes in the 

early universe, cosmological structure formation, black-hole accretion, neutron-star binaries, astrophysical 
jets, and multiphase galactic dynamics. Such an ambitious range of topics requires that the code be able 
to accurately model a wide variety of physical processes, in both the relativistic and Newtonian regimes. 
Our first contributions therefore present the results of tests designed to examine the abilities of Cosmos to 



> 



compute physical processes relevant for astrophysical problems. 

Cosmos is a collection of massively parallel, multi-dimensional, multi-physics solvers utilizing the MPI 
paradigm for parallel computing of both Newtonian and general relativistic astrophysical problems. It cur- 
rently includes several different options for computational fluid dynamics (CFD) methods, equilibrium and 
non-equilibrium primordial chemistry with 27 atomic and molecular reactions, various radiative cooling pro- 
cesses, nonequilibrium radiation flux-limited diffusion, radiation pressure, relativistic scalar fields, Newtonian 
external and self-gravity, arbitrary spacetime curvature in the form of a generic background metric, and vis- 
cous stress in a fully covariant formulation. In general, Cosmos assumes no particular symmetry in the 
equations, and is therefore designed to run on structured Cartesian meshes. However, due to the covariant 
formulation adopted for the hydrodynamics equations, Cosmos can also be run on various grid geometries 
(e.g., Cartesian, cylindrical, spherical) for problems using the relativistic fluid dynamics solvers. 

Numerical methods used to solve the hydrodynamics equations include a total variation diminishing 
(TVD) Godunov solver for Newtonian flows using Roe's (Roe 1981) approximate Riemann solver and a 
third order Runge-Kutta time-marching scheme (Shu & Osher 1988). For either Newtonian or general 
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relativistic systems, two artificial viscosity (vonNcumann & Richtmyer 1950) algorithms are available: a 
non-staggered grid method in which all variables are located at the zone-center, and a staggered grid and field 
centering method similar to the Zeus code (Stone & Norman 1992) in which scalar (vector) quantities are 
located at the zone (face) centers. Two additional options for CFD methods are included for both Newtonian 
and relativistic problems that are based on non-oscillatory central difference schemes (Jiang ct al. 1998), 
also differentiated by grid centering: staggered versus centered in time as well as space. The relativistic 
hydrodynamics algorithms have been presented in Anninos & Fragile (2003) along with several tests of 
the code, including relativistic shock tube, wall shock, and dust accretion problems. Here we emphasize 
the Newtonian and multi-physics descriptions and code tests of primordial chemistry, radiative cooling, and 
radiation diffusion coupled together with hydrodynamics. 

We proceed in §2 by describing the basic equations used by Cosmos, and in §3 by presenting tests of the 
code. The tests are designed to examine the ability of the code to follow shocks, blast waves and dynamical 
instabilities, to advect materials, to resolve heating and cooling flows, to transfer radiation, to simulate 
chemical networks, and to calculate self-gravitating gas distributions. 



2. Basic Equations 

The Newtonian multi-species mass, momentum and energy (hydrodynamic and radiation diffusion) 
continuity equations are written in an Eulerian frame as 

^+V i (pv i ) = 0, (2-1) 



dt 
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^ + ViiEv 1 ) = V, ^^—V l E^j - |vy - cpa a E + cpa p a r T\ (2-5) 

where v k is the fluid velocity assumed to be the same for each of the chemical species, e is the fluid energy 
density, p^ are the species densities satisfying p = ^ m p[ m l for the total density, <fi is the gravitational 
potential obtained from Poisson's equation V 2 = 4irGp, P is the fluid pressure, Pr — E/3 is the radiation 
pressure, a r (— 4a/ c) is the radiation constant, a is Stephan-Boltzmann's constant, c is the speed of light, E is 
the radiation energy, oy, a a , and a p are the Rosseland, absorption and Planck mean opacities, kij(T) are the 
rate coefficients for the 2-body reactions which are functions of fluid temperature T, and Ij are frequency- 
integrated photoionization and dissociation rates. The summations in (2-2) are over the N s atomic and 
molecular species included in the chemistry model, up to nine: H I, H II, He I, He II, He III, e~, H~, H2, 
H^. A total of 27 chemical reactions are included in the full network, which we summarize in Table 1 for 
convenience, but refer the reader to Abel et al. (1997) and Anninos et al. (1997) for more complete 
descriptions of the reactions, rate coefficients, and numerical methods. 

Assuming that reasonable models are provided for the mean opacities of the fluid, radiation energy 
is coupled to the net fluid momentum and internal energy, accounting for non-equilibrium heating and 
cooling effects in the single temperature and gray (spectral average) approximations. However, we have not 
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currently coupled radiation in a self-consistent fashion to the chemistry solver and ionization states of the 
fluid composition. This is a much more difficult task due to grid resolution requirements, computational time 
constraints, multi-species interactions, multi-group transport, and non-LTE effects expected to be important 
in some of our applications. We will address these problems in future developments of Cosmos. 

Three radiative cooling and heating models are implemented in the optically thin limit for A(T, pM) in 
(2-4), depending upon whether the chemistry is known. First, if chemistry is not solved, the cooling function 
is set proportional to the square of the total number density with an empirical model for the ionization 
fraction that approximates linearly the equilibrium result (2-8) described below. In particular, 
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where &m{T) is the temperature-dependent cooling rate from metals (including carbon, oxygen, neon, and 
iron lines taken from Bohringer and Hensler (1989), with a metallicity dial to scale relative to solar abun- 
dance), ej(T) is the cooling rate from hydrogen and helium lines, ju is the mass fraction of metals, fj is an 
estimate of the ionization fraction defined as min(l, max(0, (T e y — T c )/3)) with T c = leV to match roughly 
the upper and lower bounds in the equilibrium model described below. Also, 
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(2-7) 



is the gas temperature in Kelvin, /i is the mean molecular weight, m p is the proton mass, k B is Boltzmann's 
constant, and the exponential is introduced to suppress cooling at low temperature (T m j„ is typically set to 
10 4 K, with width ST ~ 500 K). Ac represents Compton cooling (or heating) due to the interaction of free 
electrons with the cosmic microwave background. 

The second model applies when the chemistry is solved in equilibrium and the ionization fraction is 
determined from equating the dominant hydrogen recombination and collisional rates (Bond et al. f984; 
Anninos et al. f 994) 

/ 
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where T e y is the temperature in electron-volts, Ih = 13.6 eV is the ionization energy of hydrogen. Assuming 
the electrons and ions have the same temperature, the gas pressure and energy can be written 
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where k (— 11605 K/eV) is a numerical constant converting between electron-volts and Kelvin. The cooling 
function takes the same form as (2-6), but in this case the ionization fraction // is computed iteratively from 
the equilibrium model fi = f of equation (2-8), and the exponential cutoff is not applied since this behavior 
is contained implicitly in /. 



-4 - 



Equilibrium chemistry is solved by setting the time derivatives in equations (2-2) to zero (but only for 
the local chemistry update, transport is still accounted for with operator splitting) and solving the following 
algebraic equations, after first making a guess for the electron number density 

nm = i I j I j / > (2-12 
«i + k 2 + k 20 /n e 

nun = n H - «hi, (2-13) 

n ^ ^ _ wm , 2 
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"Hoi = T—T—j— i (2-15) 
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"Hoin = n He ii f + 7-^- ) . (2-16) 



fc 6 A: 6 n e/ 

n e = n H n + riHeii + 2n He iii, (2-17) 

where /h is the mass fraction of hydrogen, 7Jh = /hp/w p , and y = (1 — /h)/(4/h) parameterizes the relative 
helium concentration. Equations (2-12) - (2-17) are solved iteratively in the order they are written until the 
electron density converges to a specified tolerance, typically less than 1CP 10 . The subscripts used for the 
rate coefficients ki refer to a particular chemical reaction ordered as in Table 1. 

The third cooling model is used when chemistry is solved dynamically with the full noncquilibrium 
equations (2-2). In this case, the various ionization states and concentration densities of each element 
are calculated from the chemistry equations with a stable ordered backwards differencing scheme (Anninos 
et al. 1997) and used explicitly in the cooling function as 

A(T,pM) =£ge ij (T)p[ < V bl + £ J <0V i] +tM{T){f M pf + A c , (2-18) 
»=i j=i »=i 

where ejj(T) are the cooling rates from 2-body interactions between species i and j, and Ji represents 
frequency-integrated photoionization and dissociation heating. The equation of state in this case is given by 

T=<^-^. (, 19 , 

We account for a total of eight different cooling and heating mechanisms: collisional-excitation, collisional- 
ionization, recombination, bremsstrahlung, metal- line cooling (dominantly carbon, oxygen, neon, and iron), 
molecular-hydrogen cooling, Compton cooling or heating, and photoionization heating. 

Finally we note that the following conservation equations for the chemical concentrations of hydrogen, 
helium, and charge 

Ph + Ph+ + Pa- + P H + + PH 2 = P /h, 

2 

PHc + PHc+ + PHc++ = P (1 - /h), 

1 1 1 

Ph+ - Ph- + + 4^He+ + 2^He++ = m p n e , 



are enforced after each computational cycle update. 
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3. Numerical Methods and Tests 

Cosmos uses structured, block-ordered meshes for both spatial finite differencing and finite volume dis- 
cretization methods. Depending on the hydrodynamic algorithm, state variables are either defined all at the 
zone centers (for one of the artificial viscosity schemes and for the total energy conserving methods - TVD 
Godunov and non-oscillatory central difference schemes), or on a staggered mesh for a second artificial vis- 
cosity method in which scalar and tensor quantities are zone-centered and vector variables are face-centered. 
Periodic, reflection, constant-in-time, user-specified, outgoing, and flat (vanishing first derivative) boundary 
conditions are implemented. Both the relativistic and Newtonian hydrodynamic equations are solved using 
single or multiple step time-explicit, operator-split methods with second-order spatial finite differencing. The 
gravitational potential and nonequilibrium radiation diffusion equations (discretized implicitly without the 
transport and compressive terms) are solved using a collection of linear matrix solvers from the Hypre soft- 
ware package developed at Lawrence Livermore National Laboratory (Falgout et al. 2002). Hypre includes 
options for different conjugate gradient and multigrid methods with various preconditioners to optimize 
performance. We use a multigrid solver for the gravity and diffusion tests in sections §3.5, §3.6, and §3.7. 

Since the main emphasis in this paper is on Newtonian hydrodynamics, the code tests presented in the 
following sections are designed to verify our code only in that regime, along with multi-physics (e.g. radiation 
and chemistry) coupling. We refer the interested reader to (Anninos & Fragile 2003) for discussions of the 
relativistic tests and for more explicit details of the numerical algorithms. It is anticipated that problems 
run in the Newtonian regime using Cosmos shall include microphysics (chemistry, heating, and cooling). 
Such problems are most appropriately solved with one of the artificial viscosity methods, which are written 
in an internal energy formulation. All of the hydrodynamics tests shown here are therefore computed using 
the staggered-mesh artificial viscosity method. We have, however, confirmed that results from the other 
algorithms yield comparable accuracy for most of the tests. The exceptions being the non-oscillatory central 
difference schemes which are more diffusive, in general, but particularly so for the Rayleigh- Taylor and 
spherical polytrope tests. Accuracy comparable to the artificial viscosity approach can be achieved for those 
tests if the grid resolution is roughly doubled. 

3.1. Shock Tube 

We begin testing with one of the standard problems in fluid dynamics, the shock tube or Sod problem. 
This test is characterized initially by two different fluid states separated by a membrane. At t = the 
membrane is removed and the fluid evolves in such a way that five distinct regions appear in the flow: an 
undisturbed region at each end, separated by a rarefaction wave, a contact discontinuity, and a shock wave. 
This problem provides a good test of the shock-capturing properties of the code since it has an exact solution 
(Sod 1978) against which numerical results can be compared. 

The initial state is specified by pl = 1, Pl = 1, an d Vl = to the left of the membrane, and pn = 0.125, 
Pr = 0.1, and Vr = to the right. The fluid is assumed to be an ideal gas with T = 1.4, and the integration 
domain extends over a unit grid from x = to x = 1, with the membrane located at x = 0.5. The results 
presented here were run using scalar artificial viscosity with a quadratic viscosity coefficient k q 2 = 2.0, linear 
viscosity coefficient k qi — 0.3, and Courant factor k c = 0.6, in the notation of Anninos & Fragile (2003). 
Figure 1 shows spatial profiles of the results at time t = 0.2 for the 64-zone ID case and the 64 3 -zone 3D case 
(along the main diagonal) . Table 2 summarizes the errors in p, P, and V for different grid resolutions using 
the L-l norm (i.e., ||.E(a)||i = J2i.j,k ^ x i^yj^ z k\ a i.j,k ~ ^£j,kl> wnere a i,j.k and ^i,j,k are tne numerical 
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and exact solutions, respectively, and j — k — Ayj — Azk = 1 for ID problems). The convergence rates for 
these tests are just under first order, as expected for shock-capturing methods. The slightly higher L-l norm 
errors in the 3D case are due to the fact that the error calculation is computed globally across the whole 
mesh, and so suffers from reflection effects at the grid boundaries. 

3.2. Sedov Blast Wave 

The next problem we consider is the Sedov blast wave in which energy is released at t = in the form 
of an explosion into an initially undisturbed, uniform gas. In 3D, this results in a spherical shock wave (or 
blast wave) expanding from the explosion, such that r$ oc t 2 / 5 , where r$ is the radius of the shock and t 
is the elapsed time (Sedov 1959). This problem encompasses a number of useful tests for our code, as it 
determines how well the code can follow a spherical shock wave as well as testing energy conservation. 

The initial state is specified by p = lj -Po = 3.33 x 10 -11 , and Vq = 0. The fluid is assumed to be an 
ideal gas with T — 4/3, and the integration domain extends over a unit cube. The blast wave is initiated by 
significantly increasing the energy density (relative to the background) in an approximately spherical region 
with a half- width at half-maximum of two zones and maximum cutoff radius of five zones using a Gaussian 
profile. The initial energy density contrast between the peak of the Gaussian and the background is 6.6 x 10 11 . 
The results presented here were run using a tensor artificial viscosity with a quadratic viscosity coefficient 
k q 2 = 1-0, linear viscosity coefficient k q \ — 0.3, and Courant factor k c = 0.4. Figure 2 shows the shock radius 
as a function of time fit with a curve of the form i 2 / 5 . Figure 3 shows spatial profiles of the density along 
the x, y, and z axes, as well as the diagonal at time t = 2.1 for a 64 3 -zone octant, demonstrating that the 
blast wave maintains its self-similar solution along all three major axes and the diagonal. The energy loss 
in this problem was about 19%, although much of this was at the beginning of the simulation, and the total 
energy approaches a constant value after a time t — 0.28. 

3.3. Ray leigh- Taylor Instability 

The growth of a classical Ray leigh- Taylor instability has been modeled in two-dimensions. These models 
serve as tests of the ability of the code to follow the growth of a classical instability in both the linear and 
non-linear phases, and of its ability to cleanly advect material across the grid during the non-linear growth. 

The system which we model has physical dimensions that run from to 0.1 in x, and from -0.35 to 0.15 
in z. The resolution is 128 zones in x and 1280 in z. The constant gravitational field, g = 1, points toward 
negative z. At z — 0, the density and isothermal sound speed of the dense (upper) fluid are p u = 1 and 
(c s ) 2 = 2.4, while those of the light (lower) fluid are pi = 0.1 and {c s )f = 24, such that pressure is continuous 
across the fluid interface. Away from the interface, the individual fluids are isothermal, and their densities 
vary so as to maintain hydrostatic equilibrium in the gravitational field, i.e. 

p(z) = p(O)cxp[-0(z)/c 2 ], (3-1) 

where tfi(z) is the gravitational potential, and we take <fi(0) = 0. To approximate incompressibility, the fluids 
are taken to be ideal gases, with large adiabatic index, T = 10. 

The initial distribution of the fluids, in which quantities depend only upon z, are perturbed by intro- 
ducing a vertical shift of the form 

5z{x) = Acos(27ra;/A) , (3-2) 
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whcre A is the wavelength of the perturbation, and A is its amplitude. We have examined models having 
two different values of A/X. In the first, we choose A = 0.1, and A = 0.01, such that the perturbation 
is linear, though not strongly so. Resolution requirements and computational time constraints prevent us 
from modeling a system having a substantially smaller value of A. In the second model, A = A = 0.1, such 
that this model examines the growth of initially non-linear perturbations. In introducing the perturbations, 
density and pressure are not altered, and so the fluids remain in pressure equilibrium with each other (thus 
preventing the initial growth of sound waves), but are slightly out of equilibrium with the gravitational field. 
Reflective boundary conditions are applied at the top and bottom of the grid, while periodic boundaries are 
used at the left and right. Artificial viscosity is not used in these calculations. 

Figure 4 shows snapshots of the evolution of the model having A/X — 0.1 at the times t = 0, 0.2, 0.5, 0.7, 
and 1.1. Shown, in grayscale, is the tracer material initially placed in the dense upper fluid. The initial 
perturbation can be seen in the first panel. As can be seen in the subsequent panels, the discretization 
of the grid creates short-wavelength structure superimposed upon the long-wavelength perturbation. The 
shortest wavelengths saturate quickly, but longer wavelengths persist, and can be seen superimposed upon 
the classical single-mode rollup at t = 0.5. Although the small scale vorticity spreads the tracer material 
widely across the grid, little diffusiveness is apparent. 

The early growth of the longest wavelength modes of both non-linear and linear initial perturbations 
are shown in Figure 5. Also shown is the prediction of linear theory for incompressible, constant density 
fluids, for which the instability is predicted to grow as e wt , with 

^2 = ( Pu~ Pi \ 
A \Pu + Pi) 

(Chandrasekhar 1961). The overall agreement between theory and the model with the linear initial per- 
turbation is good. The computed growth rate is found to be 14% slower than the analytic value. The slow 
growth rate may be due to the relatively large initial amplitude of the imposed perturbation, a surmise 
supported by the curvature seen in the growth rate of the numerical results. Both curves stand in contrast 
to that of the initially non-linear perturbation, which has a much smaller growth rate, as predicted by theory 
(see below). 

In the non-linear regime, the amplitude of the Rayleigh- Taylor instability is known to behave asymp- 
totically as 

A(t) = a s g^-*t* , (3-4) 

Pu + Pi 

(e.g. Youngs 1994; Glimm et al. 2001). The late-time behavior of the models is shown in Figure 6. In the 
figure are shown the penetration amplitude of the dense fluid as a function of t 2 for both the initially linear 
and initially non-linear simulations. As can be seen in the figure, both simulations display the expected 
late-time behavior, with a s ~ 0.17 for the initially linear model, and a s ~ 0.04 for the initially non-linear 
model, where a s is computed for the initial Atwood number of the fluids. The differences between the two 
models are most likely due to compressibility effects, which become more apparent at late times. 



(3-3) 
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3.4. Radiative Shock Waves 

The initial data for this test is characteristic of pre-shock flows expected from large-scale (galactic-type) 
structures 

p = 4.72 x l(T 25 g crrT 3 , 
e = 1.0 x l(T 30 g crrT 1 s~ 2 , 
v x = -u in = -1.7 x 10 7 cm s _1 , 

corresponding to a uniform flow of gas along the —x direction. Reflection boundary conditions are imposed 
at x — and we use 100 zones to resolve a spatial extent of L — 2.43 x 10~ 4 Mpc. A shock wave forms at 
the wall at x = and propagates to the right at velocity v s ~ u; n /3. As the heated gas cools radiatively, 
the shock begins to lose pressure support and slows down. Eventually the shock collapses and re-establishes 
a new pressure equilibrium closer to the wall. As gas continues to accrete, the shock front moves outwards 
again to repeat the cycle of oscillations as shown in Figure 7, where the shock position, x s , is plotted as a 
function of time in units where the grid length is set to unity and the unit of time is 10 15 seconds. 

Figure 7 shows results from two calculations: one with no chemistry in which A(T,p) oc ( o 2 T 1 / 2 , and 
a second 6-species equilibrium chemistry model that approximates the number densities of the dominant 
coolants nm, rami, «Hoi, «Hoii, ^Hciii, ne f rom equilibrium assumptions. We have also run a third calcula- 
tion to test the 6-species nonequilibrium chemistry model which explicitly solves the non-linear differential 
kinetics equations. The results in this case are nearly identical to the equilibrium calculations so we do 
not include them in Figure 7. We do, however, show the mass fraction distribution of each of the chemical 
species at the final time of the simulation for the nonequilibrium case in Figure 8. Mass fractions in the hot 
phase, where collisional ionization and recombination effects are expected to be in equilibrium, agree with 
those computed from the equilibrium model. For all these calculations we assume a perfect fluid model for 
the gas with adiabatic index r = 5/3 and a cooling function dominated by bremsstrahlung effects. 

The accuracy of these calculations is evaluated by comparing the fundamental frequency of oscillations 
to the perturbation results of Chevalier & Imamura (1982) who defined the normalized frequency as cjj = 
(2ir / P)(x s /u in ), where P is the period of the oscillations, x s is the average shock position, and u in is the 
inflow velocity. We find ujj — 0.319 and 0.315 for the no-chemistry and 6-species chemistry cases respectively. 
These results compare nicely with the perturbation estimate of 0.31. 



3.5. Marshak Wave 

In the following test, we consider the penetration of radiation from a hot source into cold material 
(Marshak 1957). Because radiative transport is very efficient in the problems considered here, significant 
penetration can occur on a timescale much shorter than the timescale for motion of the gas. We therefore 
ignore hydrodynamic transport and consider only supersonic radiation diffusion through an ambient gas in 
thermal equilibrium with the radiation and held initially at temperature To, but subject to the boundary 
condition T = T\ > Tq at one end. 

The energy diffusion equation 

^=VPV£), (3-5) 
where D = c/(3pa r ) is the diffusion coefficient, and a r — n (p / 'po) 1 (T '/T )~ m is the Rosscland mean opacity, 
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can be solved approximately for this simple case to give (Long & Tahir 1986) 



where 
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1/2 



(n + l)(n + 0.5)\ 1/2 



(3-6) 

(3-7) 
(3-8) 



n = (m + 4)/4, and ry is a numerical fitting factor of order unity. From (3-8) we clearly see that the radiation 
front location, Xf, should propagate as t 1 / 2 . 

In this work we consider two cases: a constant opacity with ay = ko, and a temperature-dependent 
opacity with ay = kq(T /To)~ m and m — 3. The computational grid is set to approximately 50 times the 
mean-free path (Z = 1/3/wy), and an appropriate stopping time tfi na i is defined using equation (3-8). The 
problem is initialized with kq = 1 cm g" - 1 , p = 1 g cm" 3 , T = 10 4 K, Ti = 10 6 K, and the grid is discretized 
with 400 zones. Figure 9 shows the profile of the radiation front for the constant opacity case normalized 
to the analytic radiation front position with r/ = 1.1. The numerical results, plotted at (0.25, 0.5, 0.75, 
l)xtfi na i, nicely display the self-similar nature of this solution. Figure 10 plots the numerical radiation 
front location as a function of time. Here we define the numerical front to be located at the half-maximum, 
although our results do not depend strongly upon this choice. The data is well fit with a t 1 / 2 curve. Figures 
11 and 12 present similar results for the temperature-dependent opacity. The analytic curve in Figure 11 is 
computed with r\ = 1.04. 



3.6. T = 2 Polytrope 

Here we test the linear matrix methods used in solving Poisson's equation for the gravitational potential 
of a compact self-gravitating source, and the ability of the code to maintain a balance between gravitational 
and pressure support forces in three dimensions. For this test we adopt an adiabatic polytrope star with 
T = 2 in spherical symmetry and hydrostatic equilibrium with solutions for the gas pressure, density and 
gravitational potential 

P = 2nGa 2 p 2 , (3-9) 
P = ^sinflY (3-10) 



r \aj 
GM 4nG Pc a 3 

Rs 



sin(l), (3-11) 



where a = (M/(47r 2 p c )) 1 / 3 for radii r < R s , and R s = an is the outer surface radius. Also, = -GM/r 
with negligible density and pressure (< 10~ 4 of the maximum central values) outside the star at radii r > R s - 
These tests are carried out for characteristic neutron star parameters with total mass M = 1.4M , central 
density p c = 2.5 x 10 15 g cm~ 3 , and radius 5.6 km. 

We ran a sequence of three simulations at different grid resolutions over two sound crossing times, 
t = 2T S = 2R S / ' yj2irGa 2 p c , with monopole boundary conditions to find mean relative errors of 0.0687, 
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0.0209, and 0.0101 for resolutions 12 3 , 24 3 , and 48 3 cells respectively. Monopole boundary conditions are 
implemented by computing the total mass Mj and mass centroid coordinates x k = J2 ceas px k AxAyAz/Mr 
in the computational domain. The potential field is then set to (f> = —GMt / \/J2k( x o ~ x c) 2 at the outer 
boundary x k — x k . Figure 13 shows spatial profiles in density for the 48 3 case along the x-axis at the initial 
time (the analytic solution), and along four separate directions at t = 2T S : the x, y, z and diagonal lines 
running through the origin. The data is displayed in dimensionless code units with length scale L = R s /N 
where N = 48 is the number of interior cells along an axis, and density scale D — M Q /L 3 . We observe no 
significant break in symmetry in the solutions, and the central peak density is maintained to good accuracy. 



3.7. r = 2 Polytrope with Radiation Pressure 

The hydrostatic polytrope solution in §3.6 can be generalized to include a radiation field and radiation 
pressure, thus providing a useful test of the coupling between gravity, fluid pressure, and radiation pressure. 
It also exercises the multiphysics operator splitting scheme in testing the ability of the code to maintain a 
balance between the three different self-consistently generated forces in three dimensions. 

The solution (3-9) - (3-11) is easily extended to include the effect of radiation pressure by assuming 
Pr = 0P, where (3 is a constant. This simplifies the solution considerably and allows for an effective 
(hydrodynamic plus radiation) pressure to counteract gravity with the same radial dependence as in the 
pure hydrodynamics case. The radiation energy equation (2-5) reduces, in the spherical, static, and thermal 
equilibrium limits, to r~ 2 d r (r 2 Dd r E) = 0, which can be solved trivially if the diffusion coefficient is set to 
D = ko/(r 2 d r E), where fco is a constant. The complete solution including radiation pressure and radially 
dependent opacity is 

p - ^ 



P 



^sinflY (3-13) 

f^sin - , 3-14 

r \aJ 



R s 

E 2irG(3a 



2 



P% (3-15) 



a r — 



3 1 + /3 

4irGc/3p c a 3 /r fr\ . fr\\ 
= fcod + fl la C ° S va)- Sm U)' (3 " 16) 

where o> = c/(3pD) is the Rosseland mean opacity, and the density and gravitational potential are un- 
changed. 

We ran a sequence of three simulations at different grid resolutions over half a sound crossing time and 
in an octant to properly specify the constant radiation boundary conditions across the star profile. The 
parameters are the same as the tests in §3.6, but with the addition of kg = 1 and f3 = 1 to give equal 
significance to the radiation and hydrodynamic pressures. We find mean relative errors of 0.148, 0.0405, and 
0.0138 for resolutions 6 3 , 12 3 , and 24 3 cells respectively, demonstrating roughly second order convergence. 
Figure 14 shows spatial profiles in density for the 24 3 case at the initial time (the analytic solution), and 
along the x, y, z and diagonal directions running through the origin at t = T s /2. 
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3.8. Astrophysical Jets 

In this last section we perform simulations of astrophysical jets as a final test of our code. A beam 
of low density material (nj = 10~ 3 cm~ 3 ) with radius Rj is injected into a homogeneous higher density 
(tia = 10nj) ambient medium along the +z axis. Within a cylindrical radius \/ x 2 + y 2 < Rj = 500 
parsecs in the z = plane, a constant inflow velocity of Vj = 10 9 cm s _1 is maintained as a boundary 
condition. Simple outflow conditions are imposed at all other external boundaries. The jet material and 
ambient medium are initialized in pressure equilibrium, with an ambient temperature of 10 7 K. An ideal 
gas equation of state with adiabatic index r = 5/3 is used for both materials. The Mach number for these 
parameters is M = Vj/c s « 36, where c s is the sound speed of the background material into which the beam 
flows. 

Figure 15 shows results at time t = 54 Rj/Vj from two calculations: one using a low resolution grid 
(64 x 64 x 128 cells) with 6 zones to cover the jet radius, and a high resolution grid (128 x 128 x 256 
cells) with 12 zones/ Rj. The physical box dimensions in both cases are set to 10Rj x 10-Rj x 20Rj. The 
results, particularly the high resolution case, clearly show all the morphological elements of astrophysical 
jets (Norman et al. 1982): a supersonic beam that ends in a bow shock, a cocoon composed of shock heated 
jet material, a working surface separating jet and shocked ambient gas, internal shock interactions, and the 
growth of Kelvin- Hclmholtz instabilities. 

Equating the ram pressure from the jet front and the equivalent pressure from the ambient medium 
yields 

V s = %=, (3-17) 

1 + VPA/PJ 

for an estimate of the velocity of the bow shock through the ambient medium, neglecting multi-dimensional 
effects. The sensitivity of the leading shock velocity V s to grid resolution is clear from Figure 15. In 
particular, we find, by tracking the bow shock position, velocities of V s = 0.32V} and 0.26Vj for the low and 
high resolution runs respectively. In comparison, equation (3-17) predicts V s = 0.24V}. Higher resolution 
allows more accurate modeling of the bow shock and working surface which effectively broadens the jet, and 
resolves to a greater extent 3D instabilities and internal shock interactions, all of which contribute to slowing 
the shock. These results are generally consistent with those of (Massaglia et al. 1996), who find jets with 
similar hypersonic Mach numbers and density ratios propagate at near (and greater than) unit efficiencies, 
defined by V s /V s > 1. 



4. Summary 

We have developed a new multidimensional, multiphysics code (Cosmos) that can be applied to a broad 
range of astrophysical problems, from highly relativistic scalar-field dominated applications in early universe 
cosmology, to black- hole accretion and multiphase Newtonian galactic dynamics. In this contribution, we 
presented the radiation-chemo- hydrodynamics equations solved in the Newtonian limit, along with numerous 
tests to gauge the accuracy and stability of the code in various multiphysics modes. In a companion paper 
(Fragile et al. 2002), Cosmos is applied to the problem of supernova-enrichment in dwarf-spheroidal galaxies, 
which utilizes many of the capabilities (robust shock capturing, radiative cooling flows, multiphase fluids, 
hydrodynamic instabilities, and gravitational potentials) presented and tested here. This work also comple- 
ments an earlier paper (Anninos & Fragile 2003) in which we presented the general relativistic equations 
solved in Cosmos. There we provided more detailed descriptions of our numerical methods for the different 
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energy and algorithmic formulations, together with numerical tests for highly relativistic hydrodynamical 
systems and black hole accretion. 

This work was performed under the auspices of the U.S. Department of Energy by University of Cali- 
fornia, Lawrence Livermore National Laboratory under Contract W-7405-Eng-48. 
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Fig. 1. — Density profiles for the 64-zone ID and 64 3 -zone 3D Sod tests covering a unit grid at time t = 0.2. 
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Fig. 2. — Shock radius as a function of time for a 64 3 -zone 3D Sedov test covering a unit grid. 
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Fig. 4. — Evolution of the Rayleigh- Taylor instability. Shown are tracer distributions of the dense upper 
fluid at times t = 0.0, 0.2, 0.5, 0.7, and 1.1. 
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Fig. 5. — Early growth rate for the Rayleigh- Taylor instability, showing the evolution of the amplitude 
normalized to its initial value versus time. Results are shown for the numerical simulations having initial 
A/X = 0.1 (linear regime; solid curve) and A/X = 1 (non-linear regime; dash-dot curve). The dashed curve 
shows the theoretical prediction for incompressible fluids in the linear regime. 











Fig. 6. — Growth rate of the Rayleigh- Taylor instability at late times. Shown is the amplitude of the dense 
fluid relative to the unperturbed interface location for both the initially linear perturbation (solid curve), 
and the initially non- linear perturbation (dot-dash curve). The amplitude is plotted as a function of t 2 , to 
demonstrate the asymptotic A(t) oc t 2 behavior. 
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Fig. 7. — Shock front position as a function of time in two radiative shock simulations: one without chemistry 
in which the cooling function is oc p 2 T 1 / 2 ; and one where the kinetics equations for a 6-species chemistry 
model with H I, H II, He I, He II, He III, and is solved in equilibrium (though we note that the 
nonequilibrium equations yield nearly identical results as the equilibrium case) . The fundamental frequency 
as defined in the text is Wj = 0.319 and 0.315 for the no-chemistry and 6-species cases respectively. Both 
compare nicely with the perturbation result of 0.31 derived by (Chevalier & Imamura 1982). 
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Fig. 8. — Mass fraction distribution (pi/p) of the species H I, H II, He I, He II, and He III at the final 
time t = 1 in the radiative shock test. The results shown are from the 6-species noncquilibrium model, 
and agree nicely with the equilibrium model in the hot post-shocked phase where collisional ionization and 
recombination balance is a good approximation. 
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Fig. 9. — Profile of the radiation front at 4 different times for a Marshak wave propagating through 
constant opacity medium. Also shown is an analytic approximation of the profile. 
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Fig. 10. — Plot of the radiation front position as a function of time for a Marshak wave propagating through 
a constant opacity medium. The data is fit with a t 1 / 2 curve, which is the expected relation. 
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Fig. 12. — Same as Figure 10 except for a temperature-dependent opacity. 
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Fig. 13. — Density profiles for a 48 3 zone test of the spherical hydrostatic r = 2 polytrope solution. Profiles 
along the x, y, z and diagonal directions through the origin are displayed at the initial time (analytic solution 
represented by the solid line), and after two sound crossing times (open circles). 
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coordinate position 

Fig. 14. — Density profiles for a 24 3 zone test of the spherical hydrostatic T = 2 polytrope solution with 
radiation pressure. Profiles along the x, y, z and diagonal directions through the origin are displayed at the 
initial time (analytic solution represented by the solid line), and after 0.5 sound crossing times (open circles). 
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Fig. 15. — Cross sections in the y — plane ol two simulations of jets propagating to the right along the 
positive z axis. From top to bottom: low resolution density, low resolution internal energy, high resolution 
density, high resolution internal energy. Six (twelve) cells are used to resolve the jet radius in the low (high) 
resolutions. 
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Table 1: Chemical gas phase reactions modeled in Cosmos, grouped by primordial versus molecular chains, 
and collisional versus photoreactive processes. The corresponding rate coefficients are referred to as h in 
the main text, where i is the reaction number defined in this table. For a more detailed description of 
the chemistry and for explicit formulas used in defining the kinetic and cooling coefficients, see Abel et al. 
(1997) and Anninos ct al. (1997). 



Table 2. L-l Norm Errors in density, pressure, and velocity for the Sod shock tube test . 



Grid TOIK |g(P)^ \\E(V)h 

32 1.62 x 10~ 2 1.67 x 10~ 2 4.82 x 10~ 2 

64 8.99 x 10~ 3 8.68 x 10~ 3 2.37 x 10~ 2 

128 4.91 x 10~ 3 4.51 x 10~ 3 1.34 x 10~ 2 

64 3 1.02 x 10~ 2 1.04 x 10~ 2 6.92 x 10~ 3 



